* Basics
clear all
set matsize 11000

capture program drop hanmerkalkan
program define hanmerkalkan, rclass

drop _all
args cor
* Set sample size
set obs 1000

* Generate e
gen e = rnormal()
gen unif = runiform()

* Generate x1
gen x1 = 2
replace x1 = 1 if unif < 0.33333333
replace x1 = 3 if unif > 0.66666666
*tab x1 //close enough
drop unif

* Generate x2 and x3
matrix M = 0, 0 // means
matrix S = 1, 1 // sd
matrix C = (1, `cor' \ `cor', 1) // correlation matrix
drawnorm x2 x3, means(M) sds(S) corr(C) 

* Generate latent-y
gen y_star = 2 + -1*x1 + 1*x2 + 0.5*x3 + e
gen y = 0
replace y = 1 if y_star > 0

* Capture mean values of each variable
qui su x1
local mean_x1 = r(mean)

qui su x2
local mean_x2 = r(mean)

qui su x3
local mean_x3 = r(mean)

* Run true probit
probit y x1 x2 x3

* Marginal effects

// AME
margins, dydx(x2 x3)
matrix define b_obs_true = r(b)
matrix define v_obs_true = r(V)
local marg_obsx2_true = b_obs_true[1,1]
return scalar marg_obsx2_true = `marg_obsx2_true'
local se_obsx2_true = sqrt(v_obs_true[1,1])
return scalar se_obsx2_true = `se_obsx2_true'
local marg_obsx3_true = b_obs_true[1,2]
return scalar marg_obsx3_true = `marg_obsx3_true'
local se_obsx3_true = sqrt(v_obs_true[2,2])
return scalar se_obsx3_true = `se_obsx3_true'

// MEM
margins, dydx(x2 x3) atmeans
matrix define b_mean_true = r(b)
matrix define v_mean_true = r(V)
local marg_meanx2_true = b_mean_true[1,1]
return scalar marg_meanx2_true = `marg_meanx2_true'
local se_meanx2_true = sqrt(v_mean_true[1,1])
return scalar se_meanx2_true = `se_meanx2_true'
local marg_meanx3_true = b_mean_true[1,2]
return scalar marg_meanx3_true = `marg_meanx3_true'
local se_meanx3_true = sqrt(v_mean_true[2,2])
return scalar se_meanx3_true = `se_meanx3_true'

* Run probit omitting x1
probit y x2 x3

* Marginal effects

// AME
margins, dydx(x2 x3)
matrix define b_obs_negx1 = r(b)
matrix define v_obs_negx1 = r(V)
local marg_obsx2_negx1 = b_obs_negx1[1,1]
return scalar marg_obsx2_negx1 = `marg_obsx2_negx1'
local se_obsx2_negx1 = sqrt(v_obs_negx1[1,1])
return scalar se_obsx2_negx1 = `se_obsx2_negx1'
local marg_obsx3_negx1 = b_obs_negx1[1,2]
return scalar marg_obsx3_negx1 = `marg_obsx3_negx1'
local se_obsx3_negx1 = sqrt(v_obs_negx1[2,2])
return scalar se_obsx3_negx1 = `se_obsx3_negx1'

// MEM
margins, dydx(x2 x3) atmeans
matrix define b_mean_negx1= r(b)
matrix define v_mean_negx1 = r(V)
local marg_meanx2_negx1 = b_mean_negx1[1,1]
return scalar marg_meanx2_negx1 = `marg_meanx2_negx1'
local se_meanx2_negx1 = sqrt(v_mean_negx1[1,1])
return scalar se_meanx2_negx1 = `se_meanx2_negx1'
local marg_meanx3_negx1 = b_mean_negx1[1,2]
return scalar marg_meanx3_negx1 = `marg_meanx3_negx1'
local se_meanx3_negx1 = sqrt(v_mean_negx1[2,2])
return scalar se_meanx3_negx1 = `se_meanx3_negx1'

* Confidence intervals
scalar lb_obsx2_negx1 = `marg_obsx2_negx1' - `se_obsx2_negx1'* invnormal(0.975)
scalar ub_obsx2_negx1 = `marg_obsx2_negx1' + `se_obsx2_negx1'* invnormal(0.975)
scalar lb_meanx2_negx1 = `marg_meanx2_negx1' - `se_meanx2_negx1'* invnormal(0.975)
scalar ub_meanx2_negx1 = `marg_meanx2_negx1' + `se_meanx2_negx1'* invnormal(0.975)
scalar lb_obsx3_negx1 = `marg_obsx3_negx1' - `se_obsx3_negx1'* invnormal(0.975)
scalar ub_obsx3_negx1 = `marg_obsx3_negx1' + `se_obsx3_negx1'* invnormal(0.975)
scalar lb_meanx3_negx1 = `marg_meanx3_negx1' - `se_meanx3_negx1'* invnormal(0.975)
scalar ub_meanx3_negx1 = `marg_meanx3_negx1' + `se_meanx3_negx1'* invnormal(0.975)

* MAB
return scalar mab_obsx2_negx1 =  abs(`marg_obsx2_negx1' - `marg_obsx2_true')
return scalar mab_obsx3_negx1 =  abs(`marg_obsx3_negx1' - `marg_obsx3_true')
return scalar mab_meanx2_negx1 = abs(`marg_meanx2_negx1' - `marg_meanx2_true')
return scalar mab_meanx3_negx1 = abs(`marg_meanx3_negx1' - `marg_meanx3_true')

* RMSE
return scalar rmse_obsx2_negx1 = (`marg_obsx2_negx1' - `marg_obsx2_true')^2
return scalar rmse_obsx3_negx1 = (`marg_obsx3_negx1' - `marg_obsx3_true')^2
return scalar rmse_meanx2_negx1 = (`marg_meanx2_negx1' - `marg_meanx2_true')^2
return scalar rmse_meanx3_negx1 = (`marg_meanx3_negx1' - `marg_meanx3_true')^2

* Coverage 
if `marg_obsx2_true' >= lb_obsx2_negx1 & `marg_obsx2_true' <= ub_obsx2_negx1 {
	return scalar cov_obsx2_negx1 = 1 
}
else {
	return scalar cov_obsx2_negx1 = 0
}

if `marg_obsx3_true' >= lb_obsx3_negx1 & `marg_obsx3_true' <= ub_obsx3_negx1 {
	return scalar cov_obsx3_negx1 = 1
}	
else {
	return scalar cov_obsx3_negx1 = 0
}	

if `marg_meanx2_true' >= lb_meanx2_negx1 & `marg_meanx2_true' <= ub_meanx2_negx1 {
	return scalar cov_meanx2_negx1 = 1 
}	
else {
 return scalar cov_meanx2_negx1 = 0
}

if `marg_meanx3_true' >= lb_meanx3_negx1 & `marg_meanx3_true' <= ub_meanx3_negx1 {
	return scalar cov_meanx3_negx1 = 1
}	
else {
	return scalar cov_meanx3_negx1 = 0
}

* Power
if 0 >= lb_obsx2_negx1 & 0 <= ub_obsx2_negx1 {
	return scalar power_obsx2_negx1 = 0 
}
else {
	return scalar power_obsx2_negx1 = 1
}

if 0 >= lb_obsx3_negx1 & 0 <= ub_obsx3_negx1 {
	return scalar power_obsx3_negx1 = 0
}	
else {
	return scalar power_obsx3_negx1 = 1
}	

if 0 >= lb_meanx2_negx1 & 0 <= ub_meanx2_negx1 {
	return scalar power_meanx2_negx1 = 0
}	
else {
 return scalar power_meanx2_negx1 = 1
}

if 0 >= lb_meanx3_negx1 & 0 <= ub_meanx3_negx1 {
	return scalar power_meanx3_negx1 = 0
}	
else {
	return scalar power_meanx3_negx1 = 1
}

* Run probit omitting x2
probit y x1 x3


* Marginal effects
* x3

// AME
margins, dydx(x3)
matrix define b_obs_negx2 = r(b)
matrix define v_obs_negx2 = r(V)
local marg_obsx3_negx2 = b_obs_negx2[1,1]
return scalar marg_obsx3_negx2 = `marg_obsx3_negx2'
local se_obsx3_negx2 = sqrt(v_obs_negx2[1,1])
return scalar se_obsx3_negx2 = `se_obsx3_negx2'

// MEM
margins, dydx(x3) atmeans
matrix define b_mean_negx2 = r(b)
matrix define v_mean_negx2 = r(V)
local marg_meanx3_negx2 = b_mean_negx2[1,1]
return scalar marg_meanx3_negx2 = `marg_meanx3_negx2'
local se_meanx3_negx2 = sqrt(v_mean_negx2[1,1])
return scalar se_meanx3_negx2 = `se_meanx3_negx2'


* Confidence intervals
scalar lb_obsx3_negx2 = `marg_obsx3_negx2' - `se_obsx3_negx2'* invnormal(0.975)
scalar ub_obsx3_negx2 = `marg_obsx3_negx2' + `se_obsx3_negx2'* invnormal(0.975)
scalar lb_meanx3_negx2 = `marg_meanx3_negx2' - `se_meanx3_negx2'* invnormal(0.975)
scalar ub_meanx3_negx2 = `marg_meanx3_negx2' + `se_meanx3_negx2'* invnormal(0.975)

* MAB
return scalar mab_obsx3_negx2 =  abs(`marg_obsx3_negx2' - `marg_obsx3_true')
return scalar mab_meanx3_negx2 = abs(`marg_meanx3_negx2' - `marg_meanx3_true')

* RMSE
return scalar rmse_obsx3_negx2 = (`marg_obsx3_negx2' - `marg_obsx3_true')^2
return scalar rmse_meanx3_negx2 = (`marg_meanx3_negx2' - `marg_meanx3_true')^2

* Coverage 

if `marg_obsx3_true' >= lb_obsx3_negx2 & `marg_obsx3_true' <= ub_obsx3_negx2 {
	return scalar cov_obsx3_negx2 = 1
}	
else {
	return scalar cov_obsx3_negx2 = 0
}	


if `marg_meanx3_true' >= lb_meanx3_negx2 & `marg_meanx3_true' <= ub_meanx3_negx2 {
	return scalar cov_meanx3_negx2 = 1
}	
else {
	return scalar cov_meanx3_negx2 = 0
}

* Power

if 0 >= lb_obsx3_negx2 & 0 <= ub_obsx3_negx2 {
	return scalar power_obsx3_negx2 = 0
}	
else {
	return scalar power_obsx3_negx2 = 1
}	

if 0 >= lb_meanx3_negx2 & 0 <= ub_meanx3_negx2 {
	return scalar power_meanx3_negx2 = 0
}	
else {
	return scalar power_meanx3_negx2 = 1
}

* Run probit omitting x3
probit y x1 x2

* Marginal effects
* x2

// AME
margins, dydx(x2)
matrix define b_obs_negx3 = r(b)
matrix define v_obs_negx3 = r(V)
local marg_obsx2_negx3 = b_obs_negx3[1,1]
return scalar marg_obsx2_negx3 = `marg_obsx2_negx3'
local se_obsx2_negx3 = sqrt(v_obs_negx3[1,1])
return scalar se_obsx2_negx3 = `se_obsx2_negx3'

// MEM
margins, dydx(x2) atmeans
matrix define b_mean_negx3 = r(b)
matrix define v_mean_negx3 = r(V)
local marg_meanx2_negx3 = b_mean_negx3[1,1]
return scalar marg_meanx2_negx3 = `marg_meanx2_negx3'
local se_meanx2_negx3 = sqrt(v_mean_negx3[1,1])
return scalar se_meanx2_negx3 = `se_meanx2_negx3'

* Confidence intervals
scalar lb_obsx2_negx3 = `marg_obsx2_negx3' - `se_obsx2_negx3'* invnormal(0.975)
scalar ub_obsx2_negx3 = `marg_obsx2_negx3' + `se_obsx2_negx3'* invnormal(0.975)
scalar lb_meanx2_negx3 = `marg_meanx2_negx3' - `se_meanx2_negx3'* invnormal(0.975)
scalar ub_meanx2_negx3 = `marg_meanx2_negx3' + `se_meanx2_negx3'* invnormal(0.975)

* MAB
return scalar mab_obsx2_negx3 =  abs(`marg_obsx2_negx3' - `marg_obsx2_true')
return scalar mab_meanx2_negx3 = abs(`marg_meanx2_negx3' - `marg_meanx2_true')

* RMSE
return scalar rmse_obsx2_negx3 = (`marg_obsx2_negx3' - `marg_obsx2_true')^2
return scalar rmse_meanx2_negx3 = (`marg_meanx2_negx3' - `marg_meanx2_true')^2

* Coverage 

if `marg_obsx2_true' >= lb_obsx2_negx3 & `marg_obsx2_true' <= ub_obsx2_negx3 {
	return scalar cov_obsx2_negx3 = 1
}	
else {
	return scalar cov_obsx2_negx3 = 0
}	


if `marg_meanx3_true' >= lb_meanx2_negx3 & `marg_meanx3_true' <= ub_meanx2_negx3 {
	return scalar cov_meanx2_negx3 = 1
}	
else {
	return scalar cov_meanx2_negx3 = 0
}

* Power

if 0 >= lb_obsx2_negx3 & 0 <= ub_obsx2_negx3 {
	return scalar power_obsx2_negx3 = 0
}	
else {
	return scalar power_obsx2_negx3 = 1
}	

if 0 >= lb_meanx2_negx3 & 0 <= ub_meanx2_negx3 {
	return scalar power_meanx2_negx3 = 0
}	
else {
	return scalar power_meanx2_negx3 = 1
}
end


